from math import *
from numpy import *
statesSize = 17
constantsSize = 36
# -----------------------------------------------------------------------
#def getInitialsAndConstants():
V_m = -70 #-82.4202; (millivolt)
C_m = 1 # (microf_per_cm2)
# i have  kept the relation between vol_myo, vol_jsr, vol_nsr and  vol_ss in
# consistent with the current paper. However, I've changed the valu of
# total volume (vol_myo). The relatoins are:
# vol_JSR/vol_myo = 0046
# vol_NSR/vol_myo = 0.0812
# vol_ss/vol_myo =  5.7469e-05
vol_myo = 4.5000e-09#25.84e-6 # (microlitre) debug
vol_JSR = 0046* vol_myo#0.12e-6# (microlitre) debug
vol_NSR = 0.0812*vol_myo# (microlitre)              debug
vol_ss =  5.7469e-05*vol_myo# (microlitre)         debug
A_cap = 5* 2.2500e-08 #1.534e-4 # (cm2) assuming the the demensions are 1.5um*1.5um and also factor 5 for effect of caviola  debug

Ca_o = 1800 # (micromolar)
R = 8.314 # (joule_per_mole_kelvin)
T = 298#  (kelvin)
F = 96.5 # (coulomb_per_millimole)

# Calmodulin constants
CMDN_tot = 50 # (micromolar)
Km_CMDN = 0.238 # (micromolar)

# Calsequestrin   constants        
# according to A.Kapla et al paper about smooth
# muscle cells and also the current paper, these numbers(CSQN_tot and Km_CSQN) about Calsequestrin are correct
CSQN_tot = 15000 # (micromolar) 
Km_CSQN = 800 # (micromolar)


# v1 is maximum RyR channel Ca2+ permeability 
v1 = 4.5        # Maximum RyR channel Ca2 permeability (per_millisecond) debug - great for controlling the pick of [Ca]i
v2 = 1.74e-5 #  (per_millisecond)
v3 = 4*0.45      #  SR Ca2-ATPase maximum pump rate (per_millisecond)  debug

tau_tr = 20    # (millisecond)
tau_xfer = 8   # (millisecond)

K_m_up = 0.5

i_CaL_max= 7  #(picoA_per_picoF)
k_plus_a = 0.006075
k_minus_a = 0.07125
k_plus_b = 0.00405
k_minus_b = 0.965
k_plus_c = 0.009
k_minus_c = 0.0008
m = 3
n = 4
E_CaL = 63
G_Ca_total = 0.1729
G_CaB = 0.000367

k_pcb = 0.0005
k_pc_max = 0.23324
k_pc_half = 20

k_NaCa = 292.8


# Initial Condition
P_O1 = 0.149102e-4     # (dimensionless)

P_O2 = 0.951726e-10   # (dimensionless)

P_C2 = 0.16774e-3       # (dimensionless)


O = 0.930308e-18        # (dimensionless)
C2 = 0.124216e-3        # (dimensionless)
C3 = 0.578679e-8        # (dimensionless)
C4 = 0.119816e-12      # (dimensionless)

I1 = 0.497923e-18       # (dimensionless)
I2 = 0.345847e-13       # (dimensionless)
I3 = 0.185106e-13       # (dimensionless)

Ca_i = 0.115001     # (micromolar)
Ca_ss = 0.115001   # (micromolar)
Ca_JSR = 1299.5     # (micromolar)
Ca_NSR = 1299.5    # (micromolar)


P_RyR = 0 #  (dimensionless)

# ICaT - T-type Ca channel
b_inf = 1.0/(1+ exp(-(V_m + 28.0)/6.1))
g_inf = 1.0/(1+ exp((V_m + 60.0)/6.6))

b = b_inf
g = g_inf     

    
initials = [ P_O1, P_O2, P_C2, O, C2, C3, C4, I1, I2, I3, Ca_i, Ca_ss, Ca_JSR, Ca_NSR, P_RyR, b_inf, g_inf]
  #  conditions = [c_m, vol_myo, vol_JSR, vol_NSR, vol_ss, A_cap, Ca_o, R, T, F, CMDN_tot, Km_CMDN, CSQN_tot, Km_CSQN, v1, v2, v3, tau_tr, tau_xfer, K_m_up, i_CaL_max, k_plus_a, k_minus_a, k_plus_b, k_minus_b, k_plus_c, k_minus_c, m, n, E_CaL, G_Ca_total, G_CaB, k_pcb, k_pc_max, k_pc_half, k_NaCa] 

   # return ( initials, conditions )
        
# -----------------------------------------------------------------------

def computeRates(t, y):
    # From: http://www.ncbi.nlm.nih.gov/pubmed/15142845
    # Title:Computer model of action potential of mouse ventricular myocytes
    # Authors:Bondarenko, Szigeti, Bett, Kim, Rasmusson, 2004
    
    global V_m, C_m, vol_myo, vol_JSR, vol_NSR ,vol_ss, A_cap, Ca_o, R, T, F, CMDN_tot, CSQN_tot, Km_CMDN, Km_CSQN, v1, v2, v3, tau_tr, tau_xfer, K_m_up,  i_CaL_max, k_plus_a, k_minus_a, k_plus_b, k_minus_b, k_plus_c, k_minus_c, m, n, E_CaL, G_CaL, G_CaB, k_pcb, k_pc_max, k_pc_half, G_Ca_total
    
    
    
    
    if t > 1500:
        V_m  = 20
    else: 
        V_m = -70 # debug
    
#    print V_m
    
    
    O = y[0]
    C2= y[1]
    C3= y[2]
    C4= y[3]
    
    P_C2= y[4]
    P_O1= y[5]
    P_O2= y[6]
    I1= y[7]
    I2= y[8] 
    I3= y[9]
    Ca_JSR= y[10]
    Ca_NSR= y[11]
    Ca_ss= y[12]
    Ca_i= y[13]
    P_RyR = y[14]
    b = y[15]
    g = y[16]
    
    
    g_max_Ca =G_Ca_total
    T_type_roleFraction = 1 # should be between 0 and 1 - zero means T-type Ca channel has no role i.e. doesn't exist
    G_CaL = (1 - T_type_roleFraction) * g_max_Ca
    G_CaT =   T_type_roleFraction * g_max_Ca # debug
    
    # Component: Ryanodine Receptor
    P_C1 = 1.00000 - (P_C2+P_O1+P_O2)
    P_O1_prime = ( k_plus_a*(Ca_ss ** n)*P_C1+ k_minus_b*P_O2+ k_minus_c*P_C2) - ( k_minus_a*P_O1+ k_plus_b*(Ca_ss ** m)*P_O1+ k_plus_c*P_O1)
    P_O2_prime =   k_plus_b*(Ca_ss ** m)*P_O1 -  k_minus_b*P_O2
    P_C2_prime =   k_plus_c*P_O1 -  k_minus_c*P_C2
    
    
    # Component: L-type Calcium Current
    
    alpha = ( 0.400000*(exp(((V_m+12.0000)/10.0000)))*((1.00000+ 0.700000*(exp(( - ((V_m+40.0000) ** 2.00000)/10.0000)))) -  0.750000*(exp(( - ((V_m+20.0000) ** 2.00000)/400.000)))))/(1.00000+ 0.120000*(exp(((V_m+12.0000)/10.0000))))
    beta =  0.0500000*(exp(( - (V_m+12.0000)/13.0000)))
    gamma = ( k_pc_max*Ca_ss)/(k_pc_half+Ca_ss)
    K_pcf =  13.0000*(1.00000 - (exp(( - ((V_m+14.5000) ** 2.00000)/100.000))))
    
    C1 = 1.00000 - (O+C2+C3+C4+I1+I2+I3)
    C2_prime = ( 4.00000*alpha*C1+ 2.00000*beta*C3) - ( beta*C2+ 3.00000*alpha*C2)
    C3_prime = ( 3.00000*alpha*C2+ 3.00000*beta*C4) - ( 2.00000*beta*C3+ 2.00000*alpha*C3)
    I1_prime = ( gamma*O+ 0.00100000*( alpha*I3 -  K_pcf*I1)+ 0.0100000*( alpha*gamma*C4 -  4.00000*beta*K_pcf*I1)) -  k_pcb*I1
    I2_prime = ( 0.00100000*( K_pcf*O -  alpha*I2)+ k_pcb*I3+ 0.00200000*( K_pcf*C4 -  4.00000*beta*I2)) -  gamma*I2
    I3_prime = ( 0.00100000*( K_pcf*I1 -  alpha*I3)+ gamma*I2+ 1.00000*gamma*K_pcf*C4) - ( 4.00000*beta*k_pcb*I3+ k_pcb*I3)
    O_prime = ( alpha*C4+ k_pcb*I1+ 0.00100000*( alpha*I2 -  K_pcf*O)) - ( 4.00000*beta*O+ gamma*O)
    C4_prime = ( 2.00000*alpha*C3+ 4.00000*beta*O+ 0.0100000*( 4.00000*k_pcb*beta*I1 -  alpha*gamma*C4)+ 0.00200000*( 4.00000*beta*I2 -  K_pcf*C4)+ 4.00000*beta*k_pcb*I3) - ( 3.00000*beta*C4+ alpha*C4+ 1.00000*gamma*K_pcf*C4)
    
    
    I_CaL =  G_CaL*O*(V_m - E_CaL)
    B_ss = (1.00000+( CMDN_tot*Km_CMDN)/((Km_CMDN+Ca_ss) ** 2.00000)) ** ( - 1.00000)
    
    # Background Ca current
    E_CaN =  (( R*T)/( 2.00000*F))*(log((Ca_o/Ca_i)))
    i_CaB =  G_CaB*(V_m - E_CaN) # (picoA_per_picoF)
    
    # Voltage dependent T-type calcium current I_CaT
    
    g_inf = 1.0/(1+ exp((V_m + 60.0)/6.6))
    alpha_g = 0.015*exp(-(V_m+71.7)/83.33)
    beta_g = 0.015*exp((V_m+71.7)/15.38)
    tau_g = 1.0/(alpha_g+ beta_g)
    
    g_prime = (g_inf -g) / tau_g
    alpha_b = 1.068*exp((V_m + 16.3)/30.0)
    beta_b =  1.068*exp(-(V_m + 16.3)/30.0)
    
    b_inf = 1.0/(1+ exp(-(V_m + 28.0)/6.1))
    tau_b = 1.0/(alpha_b+ beta_b)
    
    b_prime = (b_inf - b )/tau_b
    
    
    I_CaT =G_CaT *b*g*(V_m - 50) 
    
    I_Ca_channels= I_CaT + I_CaL
    # Component: Calcium Fluxes
    
    J_rel =  v1*(P_O1+P_O2)*(Ca_JSR - Ca_ss)*P_RyR # (micromolar_per_millisecond)
    J_tr = (Ca_NSR - Ca_JSR)/tau_tr #  (micromolar_per_millisecond)
    J_leak =  v2*(Ca_NSR - Ca_i) #  (micromolar_per_millisecond)
    J_up = ( v3*(Ca_i ** 2.00000))/((K_m_up ** 2.00000)+(Ca_i ** 2.00000)) #  (micromolar_per_millisecond)
    J_xfer = (Ca_ss - Ca_i)/tau_xfer #  (micromolar_per_millisecond)
    P_RyR_prime =  - 0.0400000*P_RyR -  (( 0.100000*I_Ca_channels)/i_CaL_max)*(exp(( - ((V_m - 5.00000) ** 2.00000)/648.000))) #  (micromolar_per_millisecond)
    
    
    
    
    # Component: Calcium Concentration
    B_i = (1.00000+( CMDN_tot*Km_CMDN)/((Km_CMDN+Ca_i) ** 2.00000)) ** ( - 1.00000)
    Ca_i_prime =  B_i*((J_leak+J_xfer) - (J_up+( ((i_CaB) )*A_cap*C_m)/( 2.00000*vol_myo*F))) # I deleted the 'J_trpn', NCX and PMCS currents terms as we are not interested in it yet
    
    Ca_ss_prime =  B_ss*(( J_rel*vol_JSR)/vol_ss - (( J_xfer*vol_myo)/vol_ss+( I_Ca_channels*A_cap*C_m)/( 2.00000*vol_ss*F)))
    
    B_JSR = (1.00000+( CSQN_tot*Km_CSQN)/((Km_CSQN+Ca_JSR) ** 2.00000)) ** ( - 1.00000)
    Ca_JSR_prime =  B_JSR*(J_tr - J_rel)
    Ca_NSR_prime = ( (J_up - J_leak)*vol_myo)/vol_NSR - ( J_tr*vol_JSR)/vol_NSR
    
    
    
    dy = zeros(17,float)

    dy[0] = O_prime
    dy[1] = C2_prime
    dy[2] = C3_prime
    dy[3] = C4_prime
    dy[4] = P_C2_prime
    dy[5] = P_O1_prime
    dy[6] = P_O2_prime
    dy[7] = I1_prime
    dy[8] = I2_prime
    dy[9] = I3_prime
    dy[10] = Ca_JSR_prime
    dy[11] = Ca_NSR_prime
    dy[12] = Ca_ss_prime
    dy[13] = Ca_i_prime
    dy[14] = P_RyR_prime
    dy[15] = b_prime
    dy[16] = g_prime
    #print y
    return  dy

# -----------------------------------------------------------------------
    
def solve_model():
    """Solve model with ODE solver"""
    from scipy.integrate import ode
    global initials 
   # Initialise constants and state variables
   # (initials, constants) = getInitialsAndConstants()

    # Set timespan to solve over
    t = linspace(0, 10, 500)
    #print 'LLLLLLLLLLLLLLLLLLLL'

    # Construct ODE object to solve
    r = ode(computeRates)
    
    r.set_integrator('vode',method='bdf',atol=1e-06, rtol=1e-06, max_step=1,with_jacobian=False)
    r.set_initial_value(initials, t[0])
   #r.set_f_params([9999])
    #print t[0]
    # Solve model
    #print 'LLLLLLLLLLLLLL'
    states = array([[0.0] * len(t)] * statesSize)
    states[:,0] = initials
    
    tEnd = 1
    dt = 1e-8
    while r.successful() and r.t < tEnd:
        #print ">>>>>>>>>>>>>>>>>>>"
        r.integrate(r.t+dt)
    
    print r.t, r.y
  

    #for time in t[1:]:
    #    if r.successful():
    #        
    #        print '>>>>>>>>>>>>>>>>>>>>', time
    #        r.integrate(time)
    #        print "*************************"
    #        states[:,time+1] = r.y
    #    else:
    #        break

    # Compute algebraic variables
    # algebraic = computeAlgebraic(constants, states, t)
    
    return (t, states)
# -----------------------------------------------------------------------
def plot_model(t, states, algebraic):
    """Plot variables against variable of integration"""
    import pylab
#    (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()
    pylab.figure(1)
    pylab.plot(t,vstack((states,algebraic)).T)
    pylab.xlabel(legend_voi)
    pylab.legend(legend_states + legend_algebraic, loc='best')
    pylab.show()

#-------------------------Main Program ----------------------------------
if __name__ == "__main__":
    (t, states) = solve_model()
 #   plot_model(t, states, algebraic)

